Shear localization in a model glass 
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Using molecular dynamics simulations, we show that a simple model of a glassy material exhibits the shear 
localization phenomenon observed in many complex fluids. At low shear rates, the system separates into a 
fluidized shear-band and an unsheared part. The two bands are characterized by a very different dynamics 
probed by a local intermediate scattering function. Furthermore, a stick-slip motion is observed at very small 
shear rates. Our results, which open the possibility of exploring complex rheological behavior using simulations, 
are compared to recent experiments on various soft glasses. 
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Shear localization is a commonly observed phenomenon in 
the rheology of complex fluids. Over some range of shear 
rates, a fluid undergoing simple shear flow in, say, a Couette 
cell tends to separate into bands parallel to the flow direction, 
with high shear rate regions coexisting with smaller shear rate 
regions. In some cases [|J, this shear-banding phenomenon 
can be understood in terms of underlying structural changes in 
the fluid, analogous to a first order phase transition. In other 
systems, however, no such changes are evident, and coexis- 
tence appears between a completely steady region (zero shear 
rate) and a sheared, fluid region. This second type of behavior 
has been observed [^, in particular in systems of 

the so-called 'soft glass' type Such systems include dense 
colloidal pastes, granular materials, emulsions, etc. Their rhe- 
ological behavior is essentially determined by the competi- 
tion between an intrinsic slow dynamics and the acceleration 
caused by the external flow |Q, g|. The large time scales 
inherent to the glassy state manifest themselves in the non- 
linear character of the rheological properties as a function of 
the shear rate 7: existence of a yield stress (the system does no 
flow until the stress a exceeds a threshold value <j y ) and non- 
linear flow curves a = (7(7), leading to a shear rate dependent 
viscosity. 

The observation of strong heterogeneities in the flow of 
such systems suggests that a global flow curve is not sufficient 
to fully characterize the flow behavior. Indeed, a simple shear- 
thinning behavior would in general imply homogeneous shear 
flow in a planar Couette cell, since the shear stress is constant 
across the cell. More generally, it remains to clarify if these 
observations in systems with very different microscopic inter- 
actions are intrinsic to glassy dynamics, that is, if a generic 
scenario for inhomogeneous shear flow can be proposed, as 
attempted in several recent studies [|l0| |TT| |. 

In order to investigate this issue, we performed molecular 
dynamics simulations of a simple and generic glass forming 
system, consisting of a binary Lennard-Jones mixture. Since 
our work is done on a simple liquid, it can be seen as a 
link between experiments on complex materials and pure phe- 
nomenology. 

The present model system has been extensively studied [|| 
Q E2p. It is known to exhibit, in the bulk state, a com- 
puter glass transition (in the sense that the relaxation time 



becomes larger than typical simulation times) at a tempera- 
ture T c ~ 0.435 for a 80:20 mixture at a density p — 1.2 (in 
Lennard-Jones units Jl2[|). All our simulation parameters are 
chosen equal to those in Refs. |^, |l2[. Since our aim is to anal- 
yse possible flow heterogeneities, we do not impose a constant 
velocity gradient over the system as was done in Ref. Q], 
where a homogeneous shear flow was imposed through the 
use of Lee-Edwards boundary conditions. Rather we confine 
the system between two solid walls, which will be driven at 
constant velocity. By doing so, we mimic an experimental 
shear cell, without imposing an homogeneous flow. 

We first equilibrate a large simulation box with periodic 
boundary conditions in all directions, at T = 0.5. The sys- 
tem is then quenched to a temperature below T c . Most of the 
results discussed in this work correspond to T = 0.2. At this 
temperature, the structural relaxation times are orders of mag- 
nitude larger than at T c . On the time scale of computer simu- 
lations, the system is in a glassy state, in which its properties 
slowly evolve with time (aging). After a time of t = 4.10 4 
[2.10 6 MD steps], we create 2 parallel solid boundaries by 
freezing all the particles outside two parallel xy-planes at po- 
sitions z wa n = ±L 2 /2. For each computer experiment, 10 
independent samples (each containing 4800 particles) are pre- 
pared using this procedure. The cell containing the fluid par- 
ticles has dimensions L x = L y = 10 and L z = 40. 

An overall shear rate is imposed by moving in the x- 
direction all the atoms of, say, the left wall (z wa u = —20) with 
a strictly constant velocity of U wa \\. This defines the total shear 
rate 7 to t = U m \\/ L z . Velocity profiles are recorded discarding 
transients (of typical duration 1/7) associated with the start of 
the shear motion. In order to remove the viscous heat created 
by the shear, the system is coupled to a heat bath. For this 
purpose, we divide the system into parallel layers of thickness 
dz = 0.25 and rescale (once every 10 integration steps) the 
y-component of the particle velocities within the layer, so as 
to impose the desired temperature T. Such a local treatment 
is necessary to keep a homogeneous temperature profile when 
flow profiles are heterogeneous. To check for a possible in- 
fluence of the thermostat, we compared, for an extremely low 
shear rate (7 tot = (7 wa u/L z = 10~ 5 ), these results with the 
output of a simulation where the inner part of the system was 
unperturbed and the walls were thermostatted instead. Both 
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methods give identical results, indicating that the results are 
not affected by the thermostat. 

Our main observation is reported in Fig. |l[ We find that, 
at low overall shear rates 7 tot , the system separates into a ho- 
mogeneously sheared band and a part which is essentially un- 
sheared. Note that, the velocity profiles are not symmetric 
with respect to the midplane. Rather, the shear band is lo- 
calized close to one of the walls, so that a single interface 
between sheared and steady region is formed. The symme- 
try between the two walls, which is a result of Galilean in- 
variance, is restored only on average, the shear band occur- 
ring equally likely on both sides of the simulation cell (see 
Fig. [l]). A similar behavior has been observed in various ex- 
periments [§,§ 1| . 

In some cases, we observed oscillations of the velocity pro- 
file between the two mentioned solutions. This effect is pos- 
sibly related to the finite size of the simulation box leading to 
a finite probability for the band to oscillate. In contrast, this 
probability would be zero for (very large) experimental sys- 
tems, thus stabilizing one of the solutions. A discussion of 
this interesting aspect is beyond the scope of this paper, and 
we postpone it to future work. 

As shown in Fig. [j], the thickness h of the sheared region 
depends on the wall velocity C/ wa n. For very small J7 wa n, h 
is of the order of a few atomic diameters and varies only 
slightly with the wall velocity. As a consequence, the shear 
rate inside the sheared region, of order U^w/h, decreases as 
Kvaii is reduced. As U wd n is further increased, h increases 
to reach the full slab size, h = L z , at a given wall velocity, 
U c . For C/waii > U c the velocity profile is therefore linear, as 
is oberved in simple fluids. These findings are in qualitative 
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FIG. 1: Filled symbols: rescaled velocity profiles, u(z) /£/ wa ii, from 
2 independent simulation runs. In both cases, the left wall is moved 
with a constant velocity U wM = 3.33 x 1CT 3 (7 to t = 0.83 x 10" 4 ). 
Due to Galilean invariance, the sheared region may be located at ei- 
ther the moving or immobile wall. Open symbols: rescaled veloc- 
ity profiles obtained at lower wall velocities [/ wa ii = 3.33 x 10~ 2 
and 3.33 x 10" 1 corresponding to overall shear rates of 7 tot = 
0.83 x 10~ 3 and 0.83 x 10" 2 . Note that the local shear rate of 
the sheared region is smaller at smaller [/wall. 



agreement with experimental results on clay suspensions [g] . 
The rather slow variation of h with respect to a change of j tot 
at small wall velocities should, however, be contrasted to re- 
ports in [^] . However, the small size of our system compared 
to the width of interfacial regions makes a more quantitative 
analysis rather difficult for the moment. 

At a given global shear rate 7 tot = 10~ 4 ([/ wa u = 0.004), 
we also studied the influence of temperature on shear localiza- 
tion. We find that increasing the temperature is qualitatively 
similar to increasing the global shear rate: the thickness of the 
sheared region is an increasing function of the temperature. 
For example, h ss 15 [see Fig. [lj] at T = 0.2 whereas h s» 25 
at T = 0.4 and h « L z (= 40) at T = 0.5. 

We find no obvious structural differences between the two 
regions of the sample. As shown in Fig. || (lower panel), static 
properties like the density profile, shear stress, and normal 
pressure are found to be constant across the film jl3[]. 

The distinction between the two bands is therefore 
purely dynamical. This can be seen for instance on the 
layer-resolved intermediate scattering function, <f> g (t;z) — 
J2i (exp[ig y (y l (t) - j/i(0))] S (z t - z)). We choose q y 
7. 1 . which corresponds to the first maximum of the static 
structure factor. The upper panel of Fig. || depicts (f> q (t;z) 
obtained from a simulation with 7 tot = 0.83 x 10~ 4 . One may 
first note the existence of two limiting behaviors of the cor- 
relation function, corresponding to the sheared and jammed 
regions, with a rather rapid change from one behavior to the 
other within a few layers. While in the jammed region (close 
to the right wall in Fig. |2|) <f> q (t; z) barely relaxes, the corre- 
lation function in the sheared region (close to the left wall) 
exhibits a two-step relaxation as observed in homogeneously 
sheared systems. This observation reflects the acceleration of 
the structural relaxation due to the flow [9]. In the jammed 
region, the system behaves as a glassy solid, and 4> q (t; z) does 
not relax to zero on the simulation time scale. 

A way to quantify the relation between the structural re- 
laxation and the variation of the shear rate across the system 
is to look at the quantity </> g (ro; z), where tq is of order of 
1 /7tot- This quantity reflects the way the system has relaxed 
on the time scale imposed by the global shear rate. Results 
for (f> q (To;z) are shown in the lower panel of Fig. The 
change of the velocity profile, when going from the sheared 
towards the unsheared region, is accompanied by a sharp jump 
in 9 (ro; z) at the interface between these two regions. The 
profile of 4>q( T o', z ) is ver y similar to that of an order param- 
eter across an interface [|l4j]. Other order parameters could 
be proposed to characterize the local dynamics: for example, 
20 the local relaxation time of <p q {t;z), defined as the time after 
which (j> q (t\ z) goes below a fraction of unity or the local dif- 
fusion coefficient parallel to the walls. All these definitions 
lead to a two-phase picture, with well defined and spatially 
constant local characteristics in both phases. 

In Fig. ||, we summarize the results by plotting the o{pf) 
flow curve of the model at T = 0.2. In this figure, we indi- 
cate the flow curve obtained for a system undergoing homoge- 
neous shear flow, taken from Ref. [^]. We also show the points 
obtained for the same system driven by the boundaries, as de- 
scribed in the present paper. Finally, we indicate the static 
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FIG. 2: Upper panel: Intermediate scattering function, <j> q (t;z), 
computed within layers of thickness dz = 2. From the bottom to 
the top, z = -17, -15, • • • , 15, 17. The temperature is T = 0.2, 
and 7 to t = 0.83 x 10 -4 . The vertical dashed line marks the time 
ro = 5754 w 0.5/7tot [this choice for To is due to larger statistical 
noise for t > 5754]. Lower panel: <t> q (To; z) (connected diamonds). 
Connected circles depict the corresponding velocity profile. Con- 
nected triangles shows the local stress a(z) and connected stars the 
local density profile p(z). 



yield stress of the system, er y . 

To obtain er y , we apply a small tangential force, Ft, acting 
on the left wall (treated as a rigid object with overdamped 
dynamics) for a certain amount of time, during which the 
velocity profile, u(z), is sampled. The force on the wall 
is then slightly increased for a new measurement before go- 
ing over to the next higher value. The static yield stress is 
then defined as the smallest force (per unit area) for which 
the average streaming velocity in the left half of the system, 
Mav, left = J -20 u(z)dz /20, exeeds it m ; n , a minimum value, 
a few times larger than the typical statistical error of tt av , left- 
We empirically find that, for a measurement time of 4.10 3 
[5.10 4 MD steps], u m i„ = 4.10~ 4 is a reasonable choice. The 
result on er y is then averaged over 10 independent runs. Again, 
for each initial configuration, the system is first equilibrated 
at T = 0.5 before being quenched to a temperature below 
T c = 0.435. This ensures the equivalence of all individual 



runs (no history dependence). The system is then propagated 
with Ft = for a time < w before the very first increment 
of the tangential force. At T = 0.2 and using an increment 
of dFj = 0.02 (once in 5.10 4 MD steps) we have obtained 
cr y = 0.596 ± 0.022, 0.658 ± 0.009 and 0.652 ± 0.015 cor- 
responding to waiting times of i w = 10 3 ,4.10 3 and 4.10 4 
respectively. Thus, already at a waiting time of t w = 4.10 3 , 
aging effects on <r y are negligible at the temperature studied. 
Note that, due to aging effects, er y might also depend on the 
speed at which the threshold value of the tangential force is 
reached. Therefore, we have determined er y for a higher value 
of dF T = 0.05 in the case of T = 0.2 and t w = 4.10 4 . This 
gives cr y = 0.66±0.015 which agrees well with 0.652±0.015 
within the error bars. 

It can be seen in Fig. [3] that er y > er(7 t0t — > 0). There- 
fore, shear-banding can be expected in the region limited by 
the vertical dotted line, which corresponds to er(7 t0t ) < er y . 
Once the yield stress cr y is added to the flow curve, the shear 
rate becomes multivalued in a range of shear stress, a situation 
encountered in several complex fluids [[IJ . This is the very ori- 
gin of the shear-banding we observe. As a consequence, this 
phenomenon should be generic for many soft glassy materials. 
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FIG. 3: Diamonds: a versus 7 under homogeneous flow conditions 
at T = 0.2. Circles: a versus 7 tot in the boundary-driven shear flow. 
The vertical dashed line is an estimate of the global shear rate below 
which shear localization is expected. 

Finally, in the very low shear rate region, we observe a time 
dependence of the shear stress characteristic of stick-slip be- 
havior, as demonstrated in Fig. This is reminiscent to what 
is observed in friction simulations [JIB]]. Due to numerical lim- 
itations, the limit between continuum sliding and stick-slip 
behavior could not be precisely located. Qualitatively, this 
behavior is obtained when the thickness of the sheared layer 
h becomes of the order of a few particle diameters, which 
also corresponds to the width of the interface separating the 
sheared and jammed regions. 

Our results lead to the following picture, (i) For global 
shear rates 7 tot smaller than a critical value, j c — U c /L z , 
the system separates spatially into two regions, one with a 
finite, approximatively uniform, shear rate, the other being 
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FIG. 4: Shear stress versus time for 7,01 = 0.83 x l(T b at T = 0.2. 
The stress rises up to a value close to <r y ~ 0.65, before suddenly 
dropping to a value smaller than the one obtained in an homogeneous 
HOW (<Tx2 — 0.4). 



jammed; (ii) when the total shear rate is increased, the thick- 
ness of the sheared region increases; (iii) for j tot > j c , h 
reaches the thickness of the slab, so that the flow is homo- 
geneous with a linear velocity profile; (iv) at very small total 
shear rates, a stick-slip phenomenon is obtained; (v) the flow 
profile closely follows the local dynamics of the material, the 
strongly sheared region corresponding to a rapid structural re- 
laxation, the jammed region to a glassy solid. 

The qualitative agreement of these phenomena with exper- 
imental observations in very different systems is remarkable: 
Ref. [|3|] describes explicitly the same scenario for the flow be- 



havior of a Laponite clay suspension, including the stick-slip 
at very small shear rates. A more quantitative study of the flow 
profiles has been performed in Ref. [|| using a Magnetic Res- 
onance Imaging technique, confirming the existence of well 
defined jammed and flow regions, whose relative width de- 
pends on the global shear rate. In contrast to molecular simu- 
lations, an experimental determination of velocity profiles to- 
gether with a local probe of the dynamics is difficult, which 
makes the observation of such effects in simulations particu- 
larly promising. 

Our results do also constrain possible phenomenological 
descriptions of the flow scenario. It is indeed important to 
remark that models relating the local shear stress to the local 
shear rate, as a = a c +aj n (a being a constant and n = 1 cor- 
responding to Bingham fluids, and n > 1 to Herschel-Buckley 
fluids), are unable to account for these results. As emphasized 
above, the shear stress is constant throughout the Couette cell, 
so that any of the afore-mentioned models would predict a 
constant shear rate in the cell, in contrast to the present re- 
sults. It appears therefore necessary to include explicitly non- 
local terms in phenomenological descriptions HQ]. The ex- 
istence of a close connection between velocity profile and lo- 
cal dynamics supports some of the assumptions of recent ap- 
proaches treating the local 'fluidity' of the system as an 'order 
parameter' [pT[]. More detailed simulations, with systematic 
investigation of size effects, should allow a direct comparison 
with the predictions of such models. 
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